#include<stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdlib.h>

#define tot (double*) malloc(sizeof(double)*totalAN);
#define totrs (double*) malloc(sizeof(double)*totalAN*rsize);
#define sd sizeof(double)
#define PI 3.14159265359
double* factorListSet(){ // OK
 double* c = (double* ) malloc(sd*1326);
 c[0]=0.2820947917738781;
 c[1]=0.4886025119029199;
 c[2]= 0.3454941494713355;
 c[3]= 0.6307831305050401;
 c[4]= 0.2575161346821264;
 c[5]= 0.1287580673410632;
 c[6]= 0.7463526651802308;
 c[7]= 0.2154534560761004;
 c[8]= 0.06813236509555216;
 c[9]= 0.02781492157551894;
 c[10]= 0.8462843753216345;
 c[11]= 0.189234939151512;
 c[12]= 0.04460310290381928;
 c[13]= 0.0119206806752224;
 c[14]= 0.004214597070904597;
 c[15]= 0.9356025796273888;
 c[16]= 0.1708168792406481;
 c[17]= 0.03228135587163618;
 c[18]= 0.006589404174225528;
 c[19]= 0.001553137458524604;
 c[20]= 0.000491145188826305;
 c[21]= 1.017107236282055;
 c[22]= 0.156943053829006;
 c[23]= 0.02481487565210345;
 c[24]= 0.00413581260868391;
 c[25]= 0.0007550926197968211;
 c[26]= 0.0001609862874555169;
 c[27]= 4.647273819914057e-05;
 c[28]= 1.092548430592079;
 c[29]= 0.1459979252047546;
 c[30]= 0.01986780112537067;
 c[31]= 0.002809731380603064;
 c[32]= 0.000423582943233983;
 c[33]= 7.059715720566384e-05;
 c[34]= 1.384524162291759e-05;
 c[35]= 3.700296470718546e-06;
 c[36]= 1.16310662292032;
 c[37]= 0.1370734300516571;
 c[38]= 0.01638340851773374;
 c[39]= 0.002016658181771346;
 c[40]= 0.0002603494517664451;
 c[41]= 3.610397299549447e-05;
 c[42]= 5.5709639801457e-06;
 c[43]= 1.017114212991521e-06;
 c[44]= 2.542785532478802e-07;
 c[45]= 1.229622689841484;
 c[46]= 0.1296136120840625;
 c[47]= 0.01381685747288017;
 c[48]= 0.001507542743711658;
 c[49]= 0.0001706956026696042;
 c[50]= 2.040202677982873e-05;
 c[51]= 2.633890331571813e-06;
 c[52]= 3.801693229872348e-07;
 c[53]= 6.519850100689623e-08;
 c[54]= 1.536743406172475e-08;
 c[55]= 1.292720736456603;
 c[56]= 0.1232560860553381;
 c[57]= 0.01186032241055152;
 c[58]= 0.001163000296325078;
 c[59]= 0.0001174807708647752;
 c[60]= 1.23835605735013e-05;
 c[61]= 1.384524162291759e-06;
 c[62]= 1.678982165396686e-07;
 c[63]= 2.284805329141706e-08;
 c[64]= 3.706443675005072e-09;
 c[65]=0.0;
 c[66]= 1.352879094951503;
 c[67]= 0.1177530108203118;
 c[68]= 0.01032762224375018;
 c[69]= 0.0009200577156223507;
 c[70]= 8.398939417550579e-05;
 c[71]= 7.936251776977327e-06;
 c[72]= 7.858060197418365e-07;
 c[73]= 8.283122738184806e-08;
 c[74]= 9.501393408603022e-09;
 c[75]= 1.226624614576148e-09;
 c[76]=0.0;
 c[77]=0.0;
 c[78]= 1.410473958869391;
 c[79]= 0.112928295511954;
 c[80]= 0.00910002138117768;
 c[81]= 0.0007430136344100779;
 c[82]= 6.191780286750648e-05;
 c[83]= 5.309407793455069e-06;
 c[84]= 4.729996402327689e-07;
 c[85]= 4.430047519252571e-08;
 c[86]= 4.430047519252572e-09;
 c[87]=0.0;
 c[88]=0.0;
 c[89]=0.0;
 c[90]=0.0;
 c[91]= 1.46580753570876;
 c[92]= 0.1086528834200811;
 c[93]= 0.008098507775955371;
 c[94]= 0.0006104479921504446;
 c[95]= 4.681922374915654e-05;
 c[96]= 3.678465622546537e-06;
 c[97]= 2.983629604350772e-07;
 c[98]= 2.521627254651016e-08;
 c[99]= 2.246444105727492e-09;
 c[100]=0.0;
 c[101]=0.0;
 c[102]=0.0;
 c[103]=0.0;
 c[104]=0.0;
 c[105]= 1.519126944936625;
 c[106]= 0.1048297184969734;
 c[107]= 0.007268633177563331;
 c[108]= 0.0005089061138323581;
 c[109]= 3.616638267560285e-05;
 c[110]= 2.623785168354005e-06;
 c[111]= 1.955653998265095e-07;
 c[112]= 1.508819816495316e-08;
 c[113]= 1.215841656708486e-09;
 c[114]=0.0;
 c[115]=0.0;
 c[116]=0.0;
 c[117]=0.0;
 c[118]=0.0;
 c[119]=0.0;
 c[120]= 1.570637328578554;
 c[121]= 0.1013842036086075;
 c[122]= 0.006571761828847014;
 c[123]= 0.0004296095103073734;
 c[124]= 2.845158486524923e-05;
 c[125]= 1.918205460301229e-06;
 c[126]= 1.323687523896326e-07;
 c[127]= 9.407037610856029e-09;
 c[128]=0.0;
 c[129]=0.0;
 c[130]=0.0;
 c[131]=0.0;
 c[132]=0.0;
 c[133]=0.0;
 c[134]=0.0;
 c[135]=0.0;
 c[136]= 1.620511203607144;
 c[137]= 0.09825792441130971;
 c[138]= 0.005979786850412427;
 c[139]= 0.000366644255972276;
 c[140]= 2.273831148908952e-05;
 c[141]= 1.432378986515563e-06;
 c[142]= 9.207680859947902e-08;
 c[143]= 6.071364964309775e-09;
 c[144]=0.0;
 c[145]=0.0;
 c[146]=0.0;
 c[147]=0.0;
 c[148]=0.0;
 c[149]=0.0;
 c[150]=0.0;
 c[151]=0.0;
 c[152]=0.0;
 c[153]= 1.668895294531136;
 c[154]= 0.09540439392102681;
 c[155]= 0.005471817261467132;
 c[156]= 0.000315915516886449;
 c[157]= 1.842456709998578e-05;
 c[158]= 1.089467491695421e-06;
 c[159]= 6.557823669848324e-08;
 c[160]= 4.036061465505948e-09;
 c[161]=0.0;
 c[162]=0.0;
 c[163]=0.0;
 c[164]=0.0;
 c[165]=0.0;
 c[166]=0.0;
 c[167]=0.0;
 c[168]=0.0;
 c[169]=0.0;
 c[170]=0.0;
 c[171]= 1.715915629394424;
 c[172]= 0.09278609064695012;
 c[173]= 0.005032032280811556;
 c[174]= 0.0002745198669795266;
 c[175]= 1.511182131758837e-05;
 c[176]= 8.421488772623757e-07;
 c[177]= 4.76772916527422e-08;
 c[178]= 2.7526497169943e-09;
 c[179]=0.0;
 c[180]=0.0;
 c[181]=0.0;
 c[182]=0.0;
 c[183]=0.0;
 c[184]=0.0;
 c[185]=0.0;
 c[186]=0.0;
 c[187]=0.0;
 c[188]=0.0;
 c[189]=0.0;
 c[190]= 1.761681409986482;
 c[191]= 0.09037234949584909;
 c[192]= 0.004648252090384827;
 c[193]= 0.0002403553935931609;
 c[194]= 1.252939101513388e-05;
 c[195]= 6.603568883778713e-07;
 c[196]= 3.529756041865882e-08;
 c[197]= 1.919934179336579e-09;
 c[198]=0.0;
 c[199]=0.0;
 c[200]=0.0;
 c[201]=0.0;
 c[202]=0.0;
 c[203]=0.0;
 c[204]=0.0;
 c[205]=0.0;
 c[206]=0.0;
 c[207]=0.0;
 c[208]=0.0;
 c[209]=0.0;
 c[210]= 1.806287998460892;
 c[211]= 0.08813782947320221;
 c[212]= 0.004310962154688459;
 c[213]= 0.0002118722309989988;
 c[214]= 1.048923844743225e-05;
 c[215]= 5.244619223716127e-07;
 c[216]= 2.65571617661055e-08;
 c[217]= 1.365953008665088e-09;
 c[218]=0.0;
 c[219]=0.0;
 c[220]=0.0;
 c[221]=0.0;
 c[222]=0.0;
 c[223]=0.0;
 c[224]=0.0;
 c[225]=0.0;
 c[226]=0.0;
 c[227]=0.0;
 c[228]=0.0;
 c[229]=0.0;
 c[230]=0.0;
 c[231]= 1.84981925508298;
 c[232]= 0.0860613804472887;
 c[233]= 0.004012632553544782;
 c[234]= 0.0001879087358370432;
 c[235]= 8.858102756970991e-06;
 c[236]= 4.213369811727295e-07;
 c[237]= 2.027158495830164e-08;
 c[238]=0.0;
 c[239]=0.0;
 c[240]=0.0;
 c[241]=0.0;
 c[242]=0.0;
 c[243]=0.0;
 c[244]=0.0;
 c[245]=0.0;
 c[246]=0.0;
 c[247]=0.0;
 c[248]=0.0;
 c[249]=0.0;
 c[250]=0.0;
 c[251]=0.0;
 c[252]=0.0;
 c[253]= 1.89234939151512;
 c[254]= 0.08412519161795824;
 c[255]= 0.003747233864656032;
 c[256]= 0.0001675813929792027;
 c[257]= 7.539843326915562e-06;
 c[258]= 3.420142387309604e-07;
 c[259]= 1.567619693034596e-08;
 c[260]=0.0;
 c[261]=0.0;
 c[262]=0.0;
 c[263]=0.0;
 c[264]=0.0;
 c[265]=0.0;
 c[266]=0.0;
 c[267]=0.0;
 c[268]=0.0;
 c[269]=0.0;
 c[270]=0.0;
 c[271]=0.0;
 c[272]=0.0;
 c[273]=0.0;
 c[274]=0.0;
 c[275]=0.0;
 c[276]= 1.933944456973762;
 c[277]= 0.08231414245617984;
 c[278]= 0.003509886827571061;
 c[279]= 0.0001502092895384877;
 c[280]= 6.463978631424566e-06;
 c[281]= 2.802490157049559e-07;
 c[282]= 1.226616618822204e-08;
 c[283]=0.0;
 c[284]=0.0;
 c[285]=0.0;
 c[286]=0.0;
 c[287]=0.0;
 c[288]=0.0;
 c[289]=0.0;
 c[290]=0.0;
 c[291]=0.0;
 c[292]=0.0;
 c[293]=0.0;
 c[294]=0.0;
 c[295]=0.0;
 c[296]=0.0;
 c[297]=0.0;
 c[298]=0.0;
 c[299]=0.0;
 c[300]= 1.974663542417147;
 c[301]= 0.0806153015433116;
 c[302]= 0.003296604831679497;
 c[303]= 0.0001352613346337005;
 c[304]= 5.578083425836789e-06;
 c[305]= 2.31617307381405e-07;
 c[306]= 9.701381442470045e-09;
 c[307]=0.0;
 c[308]=0.0;
 c[309]=0.0;
 c[310]=0.0;
 c[311]=0.0;
 c[312]=0.0;
 c[313]=0.0;
 c[314]=0.0;
 c[315]=0.0;
 c[316]=0.0;
 c[317]=0.0;
 c[318]=0.0;
 c[319]=0.0;
 c[320]=0.0;
 c[321]=0.0;
 c[322]=0.0;
 c[323]=0.0;
 c[324]=0.0;
 c[325]= 2.014559765517849;
 c[326]= 0.07901753504364936;
 c[327]= 0.00310410193677834;
 c[328]= 0.0001223187519200979;
 c[329]= 4.842645762398801e-06;
 c[330]= 1.929356253363671e-07;
 c[331]= 7.748481175203037e-09;
 c[332]=0.0;
 c[333]=0.0;
 c[334]=0.0;
 c[335]=0.0;
 c[336]=0.0;
 c[337]=0.0;
 c[338]=0.0;
 c[339]=0.0;
 c[340]=0.0;
 c[341]=0.0;
 c[342]=0.0;
 c[343]=0.0;
 c[344]=0.0;
 c[345]=0.0;
 c[346]=0.0;
 c[347]=0.0;
 c[348]=0.0;
 c[349]=0.0;
 c[350]=0.0;
 c[351]= 2.053681083307539;
 c[352]= 0.07751119753621069;
 c[353]= 0.002929647892908799;
 c[354]= 0.0001110480167994608;
 c[355]= 4.227525749644037e-06;
 c[356]= 1.618803274529695e-07;
 c[357]= 6.244669204738395e-09;
 c[358]=0.0;
 c[359]=0.0;
 c[360]=0.0;
 c[361]=0.0;
 c[362]=0.0;
 c[363]=0.0;
 c[364]=0.0;
 c[365]=0.0;
 c[366]=0.0;
 c[367]=0.0;
 c[368]=0.0;
 c[369]=0.0;
 c[370]=0.0;
 c[371]=0.0;
 c[372]=0.0;
 c[373]=0.0;
 c[374]=0.0;
 c[375]=0.0;
 c[376]=0.0;
 c[377]=0.0;
 c[378]= 2.092070967971001;
 c[379]= 0.07608788547576704;
 c[380]= 0.002770957355023711;
 c[381]= 0.0001011810566154225;
 c[382]= 3.709477511865982e-06;
 c[383]= 1.367331551655302e-07;
 c[384]= 5.074643349641621e-09;
 c[385]=0.0;
 c[386]=0.0;
 c[387]=0.0;
 c[388]=0.0;
 c[389]=0.0;
 c[390]=0.0;
 c[391]=0.0;
 c[392]=0.0;
 c[393]=0.0;
 c[394]=0.0;
 c[395]=0.0;
 c[396]=0.0;
 c[397]=0.0;
 c[398]=0.0;
 c[399]=0.0;
 c[400]=0.0;
 c[401]=0.0;
 c[402]=0.0;
 c[403]=0.0;
 c[404]=0.0;
 c[405]=0.0;
 c[406]= 2.129768972944957;
 c[407]= 0.07474023881535388;
 c[408]= 0.002626104305793813;
 c[409]= 9.250057893259113e-05;
 c[410]= 3.270389331345833e-06;
 c[411]= 1.162082247007796e-07;
 c[412]= 4.155597631240415e-09;
 c[413]=0.0;
 c[414]=0.0;
 c[415]=0.0;
 c[416]=0.0;
 c[417]=0.0;
 c[418]=0.0;
 c[419]=0.0;
 c[420]=0.0;
 c[421]=0.0;
 c[422]=0.0;
 c[423]=0.0;
 c[424]=0.0;
 c[425]=0.0;
 c[426]=0.0;
 c[427]=0.0;
 c[428]=0.0;
 c[429]=0.0;
 c[430]=0.0;
 c[431]=0.0;
 c[432]=0.0;
 c[433]=0.0;
 c[434]=0.0;
 c[435]= 2.166811210329805;
 c[436]= 0.07346178004805379;
 c[437]= 0.002493455287686542;
 c[438]= 8.482907154550921e-05;
 c[439]= 2.896016752274726e-06;
 c[440]= 9.933255511282382e-08;
 c[441]= 3.4272987602418e-09;
 c[442]=0.0;
 c[443]=0.0;
 c[444]=0.0;
 c[445]=0.0;
 c[446]=0.0;
 c[447]=0.0;
 c[448]=0.0;
 c[449]=0.0;
 c[450]=0.0;
 c[451]=0.0;
 c[452]=0.0;
 c[453]=0.0;
 c[454]=0.0;
 c[455]=0.0;
 c[456]=0.0;
 c[457]=0.0;
 c[458]=0.0;
 c[459]=0.0;
 c[460]=0.0;
 c[461]=0.0;
 c[462]=0.0;
 c[463]=0.0;
 c[464]=0.0;
 c[465]= 2.203230756026887;
 c[466]= 0.07224678259981218;
 c[467]= 0.002371616827193302;
 c[468]= 7.802046556234715e-05;
 c[469]= 2.575059075948476e-06;
 c[470]= 8.536237715559136e-08;
 c[471]= 2.845412571853045e-09;
 c[472]=0.0;
 c[473]=0.0;
 c[474]=0.0;
 c[475]=0.0;
 c[476]=0.0;
 c[477]=0.0;
 c[478]=0.0;
 c[479]=0.0;
 c[480]=0.0;
 c[481]=0.0;
 c[482]=0.0;
 c[483]=0.0;
 c[484]=0.0;
 c[485]=0.0;
 c[486]=0.0;
 c[487]=0.0;
 c[488]=0.0;
 c[489]=0.0;
 c[490]=0.0;
 c[491]=0.0;
 c[492]=0.0;
 c[493]=0.0;
 c[494]=0.0;
 c[495]=0.0;
 c[496]= 2.239057995540692;
 c[497]= 0.07109016244861498;
 c[498]= 0.002259393678604874;
 c[499]= 7.195375314183696e-05;
 c[500]= 2.298478332305522e-06;
 c[501]= 7.372372688239529e-08;
 c[502]= 2.376948353998878e-09;
 c[503]=0.0;
 c[504]=0.0;
 c[505]=0.0;
 c[506]=0.0;
 c[507]=0.0;
 c[508]=0.0;
 c[509]=0.0;
 c[510]=0.0;
 c[511]=0.0;
 c[512]=0.0;
 c[513]=0.0;
 c[514]=0.0;
 c[515]=0.0;
 c[516]=0.0;
 c[517]=0.0;
 c[518]=0.0;
 c[519]=0.0;
 c[520]=0.0;
 c[521]=0.0;
 c[522]=0.0;
 c[523]=0.0;
 c[524]=0.0;
 c[525]=0.0;
 c[526]=0.0;
 c[527]=0.0;
 c[528]= 2.274320920733615;
 c[529]= 0.06998738827464737;
 c[530]= 0.002155755395657693;
 c[531]= 6.652805584976738e-05;
 c[532]= 2.058991638214058e-06;
 c[533]= 6.396975762177483e-08;
 c[534]= 1.997105578629141e-09;
 c[535]=0.0;
 c[536]=0.0;
 c[537]=0.0;
 c[538]=0.0;
 c[539]=0.0;
 c[540]=0.0;
 c[541]=0.0;
 c[542]=0.0;
 c[543]=0.0;
 c[544]=0.0;
 c[545]=0.0;
 c[546]=0.0;
 c[547]=0.0;
 c[548]=0.0;
 c[549]=0.0;
 c[550]=0.0;
 c[551]=0.0;
 c[552]=0.0;
 c[553]=0.0;
 c[554]=0.0;
 c[555]=0.0;
 c[556]=0.0;
 c[557]=0.0;
 c[558]=0.0;
 c[559]=0.0;
 c[560]=0.0;
 c[561]= 2.309045385777095;
 c[562]= 0.06893440650860676;
 c[563]= 0.002059809370664345;
 c[564]= 6.165878068929182e-05;
 c[565]= 1.850688996628934e-06;
 c[566]= 5.574971423464406e-08;
 c[567]= 1.687063099661222e-09;
 c[568]=0.0;
 c[569]=0.0;
 c[570]=0.0;
 c[571]=0.0;
 c[572]=0.0;
 c[573]=0.0;
 c[574]=0.0;
 c[575]=0.0;
 c[576]=0.0;
 c[577]=0.0;
 c[578]=0.0;
 c[579]=0.0;
 c[580]=0.0;
 c[581]=0.0;
 c[582]=0.0;
 c[583]=0.0;
 c[584]=0.0;
 c[585]=0.0;
 c[586]=0.0;
 c[587]=0.0;
 c[588]=0.0;
 c[589]=0.0;
 c[590]=0.0;
 c[591]=0.0;
 c[592]=0.0;
 c[593]=0.0;
 c[594]=0.0;
 c[595]= 2.343255328953802;
 c[596]= 0.0679275784432832;
 c[597]= 0.001970778938358922;
 c[598]= 5.727459951756891e-05;
 c[599]= 1.668742842362807e-06;
 c[600]= 4.878614306851687e-08;
 c[601]= 1.43241048209851e-09;
 c[602]=0.0;
 c[603]=0.0;
 c[604]=0.0;
 c[605]=0.0;
 c[606]=0.0;
 c[607]=0.0;
 c[608]=0.0;
 c[609]=0.0;
 c[610]=0.0;
 c[611]=0.0;
 c[612]=0.0;
 c[613]=0.0;
 c[614]=0.0;
 c[615]=0.0;
 c[616]=0.0;
 c[617]=0.0;
 c[618]=0.0;
 c[619]=0.0;
 c[620]=0.0;
 c[621]=0.0;
 c[622]=0.0;
 c[623]=0.0;
 c[624]=0.0;
 c[625]=0.0;
 c[626]=0.0;
 c[627]=0.0;
 c[628]=0.0;
 c[629]=0.0;
 c[630]= 2.376972965719696;
 c[631]= 0.06696362717842202;
 c[632]= 0.001887985476978138;
 c[633]= 5.331505738460513e-05;
 c[634]= 1.509185376267241e-06;
 c[635]= 4.285798472617944e-08;
 c[636]= 1.222022465559807e-09;
 c[637]=0.0;
 c[638]=0.0;
 c[639]=0.0;
 c[640]=0.0;
 c[641]=0.0;
 c[642]=0.0;
 c[643]=0.0;
 c[644]=0.0;
 c[645]=0.0;
 c[646]=0.0;
 c[647]=0.0;
 c[648]=0.0;
 c[649]=0.0;
 c[650]=0.0;
 c[651]=0.0;
 c[652]=0.0;
 c[653]=0.0;
 c[654]=0.0;
 c[655]=0.0;
 c[656]=0.0;
 c[657]=0.0;
 c[658]=0.0;
 c[659]=0.0;
 c[660]=0.0;
 c[661]=0.0;
 c[662]=0.0;
 c[663]=0.0;
 c[664]=0.0;
 c[665]=0.0;
 c[666]= 2.410218957450584;
 c[667]= 0.06603959263150834;
 c[668]= 0.00181083368723047;
 c[669]= 4.97286653043065e-05;
 c[670]= 1.368736464703448e-06;
 c[671]= 3.778791414787301e-08;
 c[672]= 1.047242907767924e-09;
 c[673]=0.0;
 c[674]=0.0;
 c[675]=0.0;
 c[676]=0.0;
 c[677]=0.0;
 c[678]=0.0;
 c[679]=0.0;
 c[680]=0.0;
 c[681]=0.0;
 c[682]=0.0;
 c[683]=0.0;
 c[684]=0.0;
 c[685]=0.0;
 c[686]=0.0;
 c[687]=0.0;
 c[688]=0.0;
 c[689]=0.0;
 c[690]=0.0;
 c[691]=0.0;
 c[692]=0.0;
 c[693]=0.0;
 c[694]=0.0;
 c[695]=0.0;
 c[696]=0.0;
 c[697]=0.0;
 c[698]=0.0;
 c[699]=0.0;
 c[700]=0.0;
 c[701]=0.0;
 c[702]=0.0;
 c[703]= 2.4430125595146;
 c[704]= 0.06515279320386852;
 c[705]= 0.001738799415163349;
 c[706]= 4.647136911331539e-05;
 c[707]= 1.244669595191882e-06;
 c[708]= 3.34327780809931e-08;
 c[709]=0.0;
 c[710]=0.0;
 c[711]=0.0;
 c[712]=0.0;
 c[713]=0.0;
 c[714]=0.0;
 c[715]=0.0;
 c[716]=0.0;
 c[717]=0.0;
 c[718]=0.0;
 c[719]=0.0;
 c[720]=0.0;
 c[721]=0.0;
 c[722]=0.0;
 c[723]=0.0;
 c[724]=0.0;
 c[725]=0.0;
 c[726]=0.0;
 c[727]=0.0;
 c[728]=0.0;
 c[729]=0.0;
 c[730]=0.0;
 c[731]=0.0;
 c[732]=0.0;
 c[733]=0.0;
 c[734]=0.0;
 c[735]=0.0;
 c[736]=0.0;
 c[737]=0.0;
 c[738]=0.0;
 c[739]=0.0;
 c[740]=0.0;
 c[741]= 2.475371751684577;
 c[742]= 0.06430079296874097;
 c[743]= 0.001671419524626688;
 c[744]= 4.35053124758359e-05;
 c[745]= 1.134706714967917e-06;
 c[746]= 2.967631270159956e-08;
 c[747]=0.0;
 c[748]=0.0;
 c[749]=0.0;
 c[750]=0.0;
 c[751]=0.0;
 c[752]=0.0;
 c[753]=0.0;
 c[754]=0.0;
 c[755]=0.0;
 c[756]=0.0;
 c[757]=0.0;
 c[758]=0.0;
 c[759]=0.0;
 c[760]=0.0;
 c[761]=0.0;
 c[762]=0.0;
 c[763]=0.0;
 c[764]=0.0;
 c[765]=0.0;
 c[766]=0.0;
 c[767]=0.0;
 c[768]=0.0;
 c[769]=0.0;
 c[770]=0.0;
 c[771]=0.0;
 c[772]=0.0;
 c[773]=0.0;
 c[774]=0.0;
 c[775]=0.0;
 c[776]=0.0;
 c[777]=0.0;
 c[778]=0.0;
 c[779]=0.0;
 c[780]= 2.507313353398387;
 c[781]= 0.06348137346502772;
 c[782]= 0.001608283431176744;
 c[783]= 4.079783155811154e-05;
 c[784]= 1.036935164213738e-06;
 c[785]= 2.642355470722224e-08;
 c[786]=0.0;
 c[787]=0.0;
 c[788]=0.0;
 c[789]=0.0;
 c[790]=0.0;
 c[791]=0.0;
 c[792]=0.0;
 c[793]=0.0;
 c[794]=0.0;
 c[795]=0.0;
 c[796]=0.0;
 c[797]=0.0;
 c[798]=0.0;
 c[799]=0.0;
 c[800]=0.0;
 c[801]=0.0;
 c[802]=0.0;
 c[803]=0.0;
 c[804]=0.0;
 c[805]=0.0;
 c[806]=0.0;
 c[807]=0.0;
 c[808]=0.0;
 c[809]=0.0;
 c[810]=0.0;
 c[811]=0.0;
 c[812]=0.0;
 c[813]=0.0;
 c[814]=0.0;
 c[815]=0.0;
 c[816]=0.0;
 c[817]=0.0;
 c[818]=0.0;
 c[819]=0.0;
 c[820]= 2.538853125964903;
 c[821]= 0.06269250935154659;
 c[822]= 0.001549025990553454;
 c[823]= 3.832063337304655e-05;
 c[824]= 9.497416376052678e-07;
 c[825]= 2.359652069719277e-08;
 c[826]=0.0;
 c[827]=0.0;
 c[828]=0.0;
 c[829]=0.0;
 c[830]=0.0;
 c[831]=0.0;
 c[832]=0.0;
 c[833]=0.0;
 c[834]=0.0;
 c[835]=0.0;
 c[836]=0.0;
 c[837]=0.0;
 c[838]=0.0;
 c[839]=0.0;
 c[840]=0.0;
 c[841]=0.0;
 c[842]=0.0;
 c[843]=0.0;
 c[844]=0.0;
 c[845]=0.0;
 c[846]=0.0;
 c[847]=0.0;
 c[848]=0.0;
 c[849]=0.0;
 c[850]=0.0;
 c[851]=0.0;
 c[852]=0.0;
 c[853]=0.0;
 c[854]=0.0;
 c[855]=0.0;
 c[856]=0.0;
 c[857]=0.0;
 c[858]=0.0;
 c[859]=0.0;
 c[860]=0.0;
 c[861]= 2.570005863478459;
 c[862]= 0.06193234731236989;
 c[863]= 0.001493321497577769;
 c[864]= 3.604912065120161e-05;
 c[865]= 8.71759362326571e-07;
 c[866]= 2.113084392782858e-08;
 c[867]=0.0;
 c[868]=0.0;
 c[869]=0.0;
 c[870]=0.0;
 c[871]=0.0;
 c[872]=0.0;
 c[873]=0.0;
 c[874]=0.0;
 c[875]=0.0;
 c[876]=0.0;
 c[877]=0.0;
 c[878]=0.0;
 c[879]=0.0;
 c[880]=0.0;
 c[881]=0.0;
 c[882]=0.0;
 c[883]=0.0;
 c[884]=0.0;
 c[885]=0.0;
 c[886]=0.0;
 c[887]=0.0;
 c[888]=0.0;
 c[889]=0.0;
 c[890]=0.0;
 c[891]=0.0;
 c[892]=0.0;
 c[893]=0.0;
 c[894]=0.0;
 c[895]=0.0;
 c[896]=0.0;
 c[897]=0.0;
 c[898]=0.0;
 c[899]=0.0;
 c[900]=0.0;
 c[901]=0.0;
 c[902]=0.0;
 c[903]= 2.60078547393005;
 c[904]= 0.06119918771222799;
 c[905]= 0.001440878600036247;
 c[906]= 3.396183429840699e-05;
 c[907]= 8.018256026422763e-07;
 c[908]= 1.897313915903101e-08;
 c[909]=0.0;
 c[910]=0.0;
 c[911]=0.0;
 c[912]=0.0;
 c[913]=0.0;
 c[914]=0.0;
 c[915]=0.0;
 c[916]=0.0;
 c[917]=0.0;
 c[918]=0.0;
 c[919]=0.0;
 c[920]=0.0;
 c[921]=0.0;
 c[922]=0.0;
 c[923]=0.0;
 c[924]=0.0;
 c[925]=0.0;
 c[926]=0.0;
 c[927]=0.0;
 c[928]=0.0;
 c[929]=0.0;
 c[930]=0.0;
 c[931]=0.0;
 c[932]=0.0;
 c[933]=0.0;
 c[934]=0.0;
 c[935]=0.0;
 c[936]=0.0;
 c[937]=0.0;
 c[938]=0.0;
 c[939]=0.0;
 c[940]=0.0;
 c[941]=0.0;
 c[942]=0.0;
 c[943]=0.0;
 c[944]=0.0;
 c[945]=0.0;
 c[946]= 2.631205051777122;
 c[947]= 0.06049146858800099;
 c[948]= 0.001391435970220645;
 c[949]= 3.203999074612218e-05;
 c[950]= 7.389472842063599e-07;
 c[951]= 1.707892514789033e-08;
 c[952]=0.0;
 c[953]=0.0;
 c[954]=0.0;
 c[955]=0.0;
 c[956]=0.0;
 c[957]=0.0;
 c[958]=0.0;
 c[959]=0.0;
 c[960]=0.0;
 c[961]=0.0;
 c[962]=0.0;
 c[963]=0.0;
 c[964]=0.0;
 c[965]=0.0;
 c[966]=0.0;
 c[967]=0.0;
 c[968]=0.0;
 c[969]=0.0;
 c[970]=0.0;
 c[971]=0.0;
 c[972]=0.0;
 c[973]=0.0;
 c[974]=0.0;
 c[975]=0.0;
 c[976]=0.0;
 c[977]=0.0;
 c[978]=0.0;
 c[979]=0.0;
 c[980]=0.0;
 c[981]=0.0;
 c[982]=0.0;
 c[983]=0.0;
 c[984]=0.0;
 c[985]=0.0;
 c[986]=0.0;
 c[987]=0.0;
 c[988]=0.0;
 c[989]=0.0;
 c[990]= 2.661276943046203;
 c[991]= 0.05980775163261159;
 c[992]= 0.001344758606772777;
 c[993]= 3.026709628913576e-05;
 c[994]= 6.82273041336726e-07;
 c[995]= 1.541097711967333e-08;
 c[996]=0.0;
 c[997]=0.0;
 c[998]=0.0;
 c[999]=0.0;
 c[1000]=0.0;
 c[1001]=0.0;
 c[1002]=0.0;
 c[1003]=0.0;
 c[1004]=0.0;
 c[1005]=0.0;
 c[1006]=0.0;
 c[1007]=0.0;
 c[1008]=0.0;
 c[1009]=0.0;
 c[1010]=0.0;
 c[1011]=0.0;
 c[1012]=0.0;
 c[1013]=0.0;
 c[1014]=0.0;
 c[1015]=0.0;
 c[1016]=0.0;
 c[1017]=0.0;
 c[1018]=0.0;
 c[1019]=0.0;
 c[1020]=0.0;
 c[1021]=0.0;
 c[1022]=0.0;
 c[1023]=0.0;
 c[1024]=0.0;
 c[1025]=0.0;
 c[1026]=0.0;
 c[1027]=0.0;
 c[1028]=0.0;
 c[1029]=0.0;
 c[1030]=0.0;
 c[1031]=0.0;
 c[1032]=0.0;
 c[1033]=0.0;
 c[1034]=0.0;
 c[1035]= 2.691012803886529;
 c[1036]= 0.05914670988469412;
 c[1037]= 0.0013006346632196;
 c[1038]= 2.862862420075548e-05;
 c[1039]= 6.310703741403411e-07;
 c[1040]= 1.393801289538844e-08;
 c[1041]=0.0;
 c[1042]=0.0;
 c[1043]=0.0;
 c[1044]=0.0;
 c[1045]=0.0;
 c[1046]=0.0;
 c[1047]=0.0;
 c[1048]=0.0;
 c[1049]=0.0;
 c[1050]=0.0;
 c[1051]=0.0;
 c[1052]=0.0;
 c[1053]=0.0;
 c[1054]=0.0;
 c[1055]=0.0;
 c[1056]=0.0;
 c[1057]=0.0;
 c[1058]=0.0;
 c[1059]=0.0;
 c[1060]=0.0;
 c[1061]=0.0;
 c[1062]=0.0;
 c[1063]=0.0;
 c[1064]=0.0;
 c[1065]=0.0;
 c[1066]=0.0;
 c[1067]=0.0;
 c[1068]=0.0;
 c[1069]=0.0;
 c[1070]=0.0;
 c[1071]=0.0;
 c[1072]=0.0;
 c[1073]=0.0;
 c[1074]=0.0;
 c[1075]=0.0;
 c[1076]=0.0;
 c[1077]=0.0;
 c[1078]=0.0;
 c[1079]=0.0;
 c[1080]=0.0;
 c[1081]= 2.720423653362309;
 c[1082]= 0.05850711688397001;
 c[1083]= 0.001258872718479037;
 c[1084]= 2.711174328577707e-05;
 c[1085]= 5.847068929770086e-07;
 c[1086]= 1.263363950714379e-08;
 c[1087]=0.0;
 c[1088]=0.0;
 c[1089]=0.0;
 c[1090]=0.0;
 c[1091]=0.0;
 c[1092]=0.0;
 c[1093]=0.0;
 c[1094]=0.0;
 c[1095]=0.0;
 c[1096]=0.0;
 c[1097]=0.0;
 c[1098]=0.0;
 c[1099]=0.0;
 c[1100]=0.0;
 c[1101]=0.0;
 c[1102]=0.0;
 c[1103]=0.0;
 c[1104]=0.0;
 c[1105]=0.0;
 c[1106]=0.0;
 c[1107]=0.0;
 c[1108]=0.0;
 c[1109]=0.0;
 c[1110]=0.0;
 c[1111]=0.0;
 c[1112]=0.0;
 c[1113]=0.0;
 c[1114]=0.0;
 c[1115]=0.0;
 c[1116]=0.0;
 c[1117]=0.0;
 c[1118]=0.0;
 c[1119]=0.0;
 c[1120]=0.0;
 c[1121]=0.0;
 c[1122]=0.0;
 c[1123]=0.0;
 c[1124]=0.0;
 c[1125]=0.0;
 c[1126]=0.0;
 c[1127]=0.0;
 c[1128]= 2.749519921161698;
 c[1129]= 0.05788783709042469;
 c[1130]= 0.001219299419741175;
 c[1131]= 2.570508877402524e-05;
 c[1132]= 5.426348493370343e-07;
 c[1133]= 1.147550435625091e-08;
 c[1134]=0.0;
 c[1135]=0.0;
 c[1136]=0.0;
 c[1137]=0.0;
 c[1138]=0.0;
 c[1139]=0.0;
 c[1140]=0.0;
 c[1141]=0.0;
 c[1142]=0.0;
 c[1143]=0.0;
 c[1144]=0.0;
 c[1145]=0.0;
 c[1146]=0.0;
 c[1147]=0.0;
 c[1148]=0.0;
 c[1149]=0.0;
 c[1150]=0.0;
 c[1151]=0.0;
 c[1152]=0.0;
 c[1153]=0.0;
 c[1154]=0.0;
 c[1155]=0.0;
 c[1156]=0.0;
 c[1157]=0.0;
 c[1158]=0.0;
 c[1159]=0.0;
 c[1160]=0.0;
 c[1161]=0.0;
 c[1162]=0.0;
 c[1163]=0.0;
 c[1164]=0.0;
 c[1165]=0.0;
 c[1166]=0.0;
 c[1167]=0.0;
 c[1168]=0.0;
 c[1169]=0.0;
 c[1170]=0.0;
 c[1171]=0.0;
 c[1172]=0.0;
 c[1173]=0.0;
 c[1174]=0.0;
 c[1175]=0.0;
 c[1176]= 2.778311490808207;
 c[1177]= 0.05728781739681247;
 c[1178]= 0.00118175744029733;
 c[1179]= 2.43985682194388e-05;
 c[1180]= 5.043783222424648e-07;
 c[1181]= 1.044460787186219e-08;
 c[1182]=0.0;
 c[1183]=0.0;
 c[1184]=0.0;
 c[1185]=0.0;
 c[1186]=0.0;
 c[1187]=0.0;
 c[1188]=0.0;
 c[1189]=0.0;
 c[1190]=0.0;
 c[1191]=0.0;
 c[1192]=0.0;
 c[1193]=0.0;
 c[1194]=0.0;
 c[1195]=0.0;
 c[1196]=0.0;
 c[1197]=0.0;
 c[1198]=0.0;
 c[1199]=0.0;
 c[1200]=0.0;
 c[1201]=0.0;
 c[1202]=0.0;
 c[1203]=0.0;
 c[1204]=0.0;
 c[1205]=0.0;
 c[1206]=0.0;
 c[1207]=0.0;
 c[1208]=0.0;
 c[1209]=0.0;
 c[1210]=0.0;
 c[1211]=0.0;
 c[1212]=0.0;
 c[1213]=0.0;
 c[1214]=0.0;
 c[1215]=0.0;
 c[1216]=0.0;
 c[1217]=0.0;
 c[1218]=0.0;
 c[1219]=0.0;
 c[1220]=0.0;
 c[1221]=0.0;
 c[1222]=0.0;
 c[1223]=0.0;
 c[1224]=0.0;
 c[1225]= 2.806807738882166;
 c[1226]= 0.05670607959001315;
 c[1227]= 0.00114610370472699;
 c[1228]= 2.318319646198824e-05;
 c[1229]= 4.695225605470426e-07;
 c[1230]= 9.524744358118536e-09;
 c[1231]=0.0;
 c[1232]=0.0;
 c[1233]=0.0;
 c[1234]=0.0;
 c[1235]=0.0;
 c[1236]=0.0;
 c[1237]=0.0;
 c[1238]=0.0;
 c[1239]=0.0;
 c[1240]=0.0;
 c[1241]=0.0;
 c[1242]=0.0;
 c[1243]=0.0;
 c[1244]=0.0;
 c[1245]=0.0;
 c[1246]=0.0;
 c[1247]=0.0;
 c[1248]=0.0;
 c[1249]=0.0;
 c[1250]=0.0;
 c[1251]=0.0;
 c[1252]=0.0;
 c[1253]=0.0;
 c[1254]=0.0;
 c[1255]=0.0;
 c[1256]=0.0;
 c[1257]=0.0;
 c[1258]=0.0;
 c[1259]=0.0;
 c[1260]=0.0;
 c[1261]=0.0;
 c[1262]=0.0;
 c[1263]=0.0;
 c[1264]=0.0;
 c[1265]=0.0;
 c[1266]=0.0;
 c[1267]=0.0;
 c[1268]=0.0;
 c[1269]=0.0;
 c[1270]=0.0;
 c[1271]=0.0;
 c[1272]=0.0;
 c[1273]=0.0;
 c[1274]=0.0;
 c[1275]= 2.835017570693472;
 c[1276]= 0.05614171363835682;
 c[1277]= 0.00111220784184347;
 c[1278]= 2.205095481555296e-05;
 c[1279]= 4.377050834318466e-07;
 c[1280]= 8.702045102005835e-09;
 c[1281]=0.0;
 c[1282]=0.0;
 c[1283]=0.0;
 c[1284]=0.0;
 c[1285]=0.0;
 c[1286]=0.0;
 c[1287]=0.0;
 c[1288]=0.0;
 c[1289]=0.0;
 c[1290]=0.0;
 c[1291]=0.0;
 c[1292]=0.0;
 c[1293]=0.0;
 c[1294]=0.0;
 c[1295]=0.0;
 c[1296]=0.0;
 c[1297]=0.0;
 c[1298]=0.0;
 c[1299]=0.0;
 c[1300]=0.0;
 c[1301]=0.0;
 c[1302]=0.0;
 c[1303]=0.0;
 c[1304]=0.0;
 c[1305]=0.0;
 c[1306]=0.0;
 c[1307]=0.0;
 c[1308]=0.0;
 c[1309]=0.0;
 c[1310]=0.0;
 c[1311]=0.0;
 c[1312]=0.0;
 c[1313]=0.0;
 c[1314]=0.0;
 c[1315]=0.0;
 c[1316]=0.0;
 c[1317]=0.0;
 c[1318]=0.0;
 c[1319]=0.0;
 c[1320]=0.0;
 c[1321]=0.0;
 c[1322]=0.0;
 c[1323]=0.0;
 c[1324]=0.0;
 c[1325]=0.0;
  return c;
}
double* getws(){ // OK
 double* c = (double* ) malloc(sd*100);
c[0] = 7.34634490505672E-4;
c[1] = 0.001709392653518105;
c[2] = 0.002683925371553482;
c[3] = 0.003655961201326375;
c[4] = 0.00462445006342212;
c[5] = 0.00558842800386552;
c[6] = 0.00654694845084532;
c[7] = 0.007499073255464712;
c[8] = 0.008443871469668971;
c[9] = 0.00938041965369446;
c[10] = 0.01030780257486897;
c[11] = 0.01122511402318598;
c[12] = 0.012131457662979497;
c[13] = 0.0130259478929715423;
c[14] = 0.0139077107037187727;
c[15] = 0.014775884527441302;
c[16] = 0.015629621077546003;
c[17] = 0.01646808617614521;
c[18] = 0.017290460568323582;
c[19] = 0.018095940722128117;
c[20] = 0.018883739613374905;
c[21] = 0.019653087494435306;
c[22] = 0.020403232646209433;
c[23] = 0.021133442112527642;
c[24] = 0.02184300241624739;
c[25] = 0.022531220256336273;
c[26] = 0.023197423185254122;
c[27] = 0.023840960265968206;
c[28] = 0.02446120270795705;
c[29] = 0.02505754448157959;
c[30] = 0.02562940291020812;
c[31] = 0.02617621923954568;
c[32] = 0.026697459183570963;
c[33] = 0.02719261344657688;
c[34] = 0.027661198220792388;
c[35] = 0.028102755659101173;
c[36] = 0.028516854322395098;
c[37] = 0.028903089601125203;
c[38] = 0.029261084110638277;
c[39] = 0.02959048805991264;
c[40] = 0.029890979593332831;
c[41] = 0.03016226510516914;
c[42] = 0.03040407952645482;
c[43] = 0.03061618658398045;
c[44] = 0.03079837903115259;
c[45] = 0.030950478850490988;
c[46] = 0.031072337427566517;
c[47] = 0.031163835696209907;
c[48] = 0.03122488425484936;
c[49] = 0.03125542345386336;
c[50] = 0.031255423453863357;
c[51] = 0.03122488425484936;
c[52] = 0.0311638356962099068;
c[53] = 0.031072337427566517;
c[54] = 0.030950478850490988;
c[55] = 0.03079837903115259;
c[56] = 0.030616186583980449;
c[57] = 0.03040407952645482;
c[58] = 0.0301622651051691449;
c[59] = 0.02989097959333283;
c[60] = 0.029590488059912643;
c[61] = 0.029261084110638277;
c[62] = 0.028903089601125203;
c[63] = 0.0285168543223951;
c[64] = 0.02810275565910117;
c[65] = 0.02766119822079239;
c[66] = 0.02719261344657688;
c[67] = 0.02669745918357096;
c[68] = 0.02617621923954568;
c[69] = 0.025629402910208116;
c[70] = 0.02505754448157959;
c[71] = 0.024461202707957053;
c[72] = 0.02384096026596821;
c[73] = 0.023197423185254122;
c[74] = 0.0225312202563362727;
c[75] = 0.021843002416247386;
c[76] = 0.02113344211252764;
c[77] = 0.020403232646209433;
c[78] = 0.019653087494435306;
c[79] = 0.0188837396133749046;
c[80] = 0.018095940722128117;
c[81] = 0.017290460568323582;
c[82] = 0.016468086176145213;
c[83] = 0.015629621077546003;
c[84] = 0.0147758845274413;
c[85] = 0.013907710703718773;
c[86] = 0.01302594789297154;
c[87] = 0.0121314576629795;
c[88] = 0.01122511402318598;
c[89] = 0.01030780257486897;
c[90] = 0.00938041965369446;
c[91] = 0.008443871469668971;
c[92] = 0.007499073255464712;
c[93] = 0.00654694845084532;
c[94] = 0.005588428003865515;
c[95] = 0.00462445006342212;
c[96] = 0.003655961201326375;
c[97] = 0.002683925371553482;
c[98] = 0.001709392653518105;
c[99] = 7.3463449050567E-4;
return c;
}

//=========================================================
double factorY(int l, int m, double* c){ // OK
  return  c[(l*(l+1))/2 + m];//l+1
}
//=========================================================
double* getoOr(double* r, int rsize){
  double* oOr = (double*) malloc(sd*rsize);
  for(int w = 0; w < rsize; w++){
    oOr[w] = 1/r[w];
  }
  return oOr;
}
//=========================================================
double* getrw2(double* r, int rsize){
  double* rw2 = (double*) malloc(sd*rsize);
  for(int w = 0; w < rsize; w++){
    rw2[w] = r[w]*r[w];
  }
  return rw2;
}
//=========================================================
void expMs(double* rExpDiff, double alpha, double* r, double* ri, int isize, int rsize){
  double rDiff;
  for(int i = 0; i < isize; i++){
    for(int w = 0; w < rsize; w++){
      rDiff = r[w] - ri[i];
    if(rDiff > 5.0 ){rExpDiff[rsize*i + w] = 0.0;}
    else {rExpDiff[rsize*i + w] = exp(-alpha*rDiff*rDiff);}
    }
  }
}
//=========================================================
void expPs(double* rExpSum, double alpha, double* r, double* ri, int isize, int rsize){
  double rSum;
  for(int i = 0; i < isize; i++){
    for(int w = 0; w < rsize; w++){
    rSum = r[w] + ri[i];
    if(rSum > 5.0 ){rExpSum[rsize*i + w] = 0.0;}
    else {rExpSum[rsize*i + w] = exp(-alpha*rSum*rSum);}
    }
  }
}
//=========================================================
int getFilteredPos(double* x, double* y, double* z,double* xNow, double* yNow, double* zNow, double* ri, double* rw, double rCut, double* oOri, double* oO4arri, double* minExp, double* pluExp,int* isCenter, double alpha, double* Apos, double* Hpos,int* typeNs, int rsize, int Ihpos, int Itype){//OK

  int shiftType = 0;
  int icount = 0;
  double ri2;
  double oOa = 1/alpha;
  double Xi; double Yi; double Zi;
  double* oO4ari = (double*) malloc(sd*typeNs[Itype]);

  for(int i = 0; i < Itype ; i++){
    shiftType += typeNs[i];
  }

  for(int i = 0; i < typeNs[Itype]; i++){
    Xi = Apos[3*shiftType + 3*i    ] - Hpos[3*Ihpos    ];
    Yi = Apos[3*shiftType + 3*i + 1] - Hpos[3*Ihpos + 1];
    Zi = Apos[3*shiftType + 3*i + 2] - Hpos[3*Ihpos + 2];
    ri2 = Xi*Xi + Yi*Yi + Zi*Zi;
    if(ri2<=1e-12) isCenter[0] = 1;
    if(ri2 < (rCut + 5)*(rCut + 5) && ri2 > 1e-12){ // 25 -> halo +5 Ang
      ri[icount] = sqrt(ri2);
      xNow[icount] = Xi; yNow[icount] = Yi; zNow[icount] = Zi;
      oOri[icount] = 1/ri[icount];
      oO4ari[icount] = 0.25*oOa*oOri[icount];
      icount++;
    }
  }
  //countMax = isize ----------------------------
  double* oOr = getoOr(rw, rsize);
  for(int i = 0; i < icount; i++){
    for(int w = 0; w < rsize; w++){
      oO4arri[rsize*i + w] = oO4ari[i]*oOr[w];
    }
  }
  expMs(minExp,alpha,rw,ri,icount,rsize);
  expPs(pluExp,alpha,rw,ri,icount,rsize);

//  free(ri);
  free(oO4ari);

  return icount;
}
//=========================================================
double* getFlir(double* oO4arri,double* ri, double* minExp, double* pluExp, int icount, int rsize, int lMax){//OK
  double* Flir = (double*) malloc(sd*(lMax+1)*icount*rsize);
//  double* rw =   getrw(100, 6);
//  int count = 0;
  //l=0
  for(int i = 0; i < icount; i++){
///    if(ri[i] < 0.01){}
    for(int w = 0; w < rsize; w++){
      Flir[rsize*i + w] = oO4arri[rsize*i + w]*(minExp[rsize*i + w] - pluExp[rsize*i + w]);
 //    exit(1);
    }
  }
  //l=1
  if(lMax>0)
  for(int i = 0; i < icount; i++){
///    if(ri[i] < 0.01){}
    for(int w = 0; w < rsize; w++){
      Flir[rsize*icount + rsize*i + w] = oO4arri[rsize*i + w]*(minExp[rsize*i + w] + pluExp[rsize*i + w] - 2*Flir[rsize*i + w]);
    }
  }
  //l>1
  if(lMax>1)
  for(int l = 2; l < lMax+1; l++){
    for(int i = 0; i < icount; i++){
 ///   if(ri[i] < 0.01){}
      for(int w = 0; w < rsize; w++){
        Flir[l*rsize*icount+rsize*i+w] = Flir[(l-2)*rsize*icount+rsize*i+w] - oO4arri[rsize*i+w]*(4*l-2)*Flir[(l-1)*rsize*icount+rsize*i+w] ;
   if(Flir[l*rsize*icount+rsize*i+w] < 0) Flir[l*rsize*icount+rsize*i+w]=0.0; // Very Important!!!

//    }
      }
    }
  }


  return Flir;

}
//=========================================================
double legendre_poly(int l, int m, double x){ // OK

  double fact,pll,pmm,pmmp1,somx2;
  int ll;

  if (m < 0 || m > l || fabs(x) > 1.0){ printf("ERROR: Bad arguments in routine legendre_poly"); exit(1);}

  pmm = 1.0;

  if(m > 0) {
    somx2=sqrt((1.0 - x)*(1.0 + x));
    fact=1.0;
    for(int i=1; i <= m; i++)
        {
          pmm *= -fact*somx2;
          fact += 2.0;
           }
   }

  if(l == m) return pmm;

  else{
    pmmp1 = x*(2*m+1)*pmm;
    if(l==(m+1)) return pmmp1;
    else{
      for(ll=m+2; ll<=l; ll++){
        pll=(x*(2*ll-1)*pmmp1 - (ll+m-1)*pmm)/ (double) (ll-m);
        pmm = pmmp1;
        pmmp1= pll;

      }
      return pll;
    }
  }

}
//=========================================================
double* getYlmi(double* x, double* y, double* z, double* oOri, double* cf, int icount, int lMax){ // OK
  double* Ylmi = (double*) malloc(2*sd*(lMax+1)*(lMax+1)*icount);
  double* legPol = (double*) malloc(sd*(lMax+1)*(lMax+1)*icount);
  double* ChiCos= (double*) malloc(sd*(lMax+1)*icount);
  double* ChiSin= (double*) malloc(sd*(lMax+1)*icount);
  double myAtan2;

  for(int i = 0; i < icount; i++){
    for(int l = 0; l < lMax + 1; l++){
      for(int m = 0; m < l+1; m++){
        legPol[icount*(lMax+1)*l + icount*m + i] = legendre_poly(l,m,z[i]*oOri[i]);
      }
    }

    for(int m = 0; m < lMax+1; m++){
      myAtan2 = m*atan2(y[i],x[i]);
      ChiCos[m*icount + i] = cos(myAtan2);
      ChiSin[m*icount + i] = sin(myAtan2);
//if(y<=0){ ChiSin[m*icount + i] = -sqrt(1 - ChiCos[m*icount + i]*ChiCos[m*icount + i]);}
//else{ChiSin[m*icount + i] = sqrt(1 - ChiCos[m*icount + i]*ChiCos[m*icount + i]);}
    }
  }

  for(int l = 0; l < lMax+1; l++){
    for(int m = 0; m < l+1; m++){//l+1
      for(int i = 0; i < icount; i++){

        Ylmi[2*(lMax+1)*icount*l + 2*icount*m + 2*i]
          =  factorY(l,m,cf)*legPol[icount*(lMax+1)*l + icount*m + i]*ChiCos[m*icount + i];
        Ylmi[2*(lMax+1)*icount*l + 2*icount*m + 2*i + 1]
          = factorY(l,m,cf)*legPol[icount*(lMax+1)*l + icount*m + i]*ChiSin[m*icount + i];

      }
    }
  }
  free(legPol); free(ChiCos); free(ChiSin);

  return Ylmi;
}
//=========================================================
double* getIntegrand(double* Flir, double* Ylmi,int rsize, int icount, int lMax){

  double* summed = (double*) malloc(2*sd*(lMax+1)*rsize*(lMax+1));
  double realY;
  double imagY;

  for(int i = 0; i < 2*(lMax+1)*rsize*(lMax+1); i++){summed[i] = 0.0;}

  for(int l = 0; l < lMax+1; l++){
    double summe = 0;
    for(int m = 0; m < l+1; m++){//l+1
     for(int i = 0; i < icount; i++){
      realY = Ylmi[2*(lMax+1)*icount*l + 2*icount*m + 2*i    ];
      imagY = Ylmi[2*(lMax+1)*icount*l + 2*icount*m + 2*i  + 1 ];
      for(int rw = 0; rw < rsize; rw++){
         summed[2*(lMax+1)*l*rsize + 2*m*rsize + 2*rw    ] += Flir[l*rsize*icount + rsize*i + rw] * realY;
         summed[2*(lMax+1)*l*rsize + 2*m*rsize + 2*rw + 1] += Flir[l*rsize*icount + rsize*i + rw] * imagY;
      }

    }
  }
  }
  return summed;
}
//=========================================================
void getC(double* Cs, double* ws, double* rw2, double * gns, double* summed, double rCut,int lMax, int rsize, int gnsize,int* isCenter, double alpha){

  for(int i = 0; i < 2*(lMax+1)*(lMax+1)*gnsize; i++){ Cs[i] = 0.0;}
  double  theSummedValue = 0;

  for(int n = 0; n < gnsize; n++){
    //for i0 case
    if(isCenter[0]==1){
      for(int rw = 0; rw < rsize; rw++){
        Cs[2*(lMax+1)*(lMax+1)*n] += 0.5*0.564189583547756*rw2[rw]*ws[rw]*gns[rsize*n + rw]*exp(-alpha*rw2[rw]);
      }
    }
    for(int l = 0; l < lMax+1; l++){
      for(int m = 0; m < l+1; m++){ // l+1
        for(int rw = 0; rw < rsize; rw++){

          Cs[2*(lMax+1)*(lMax+1)*n + l*2*(lMax+1) + 2*m    ] += rw2[rw]*ws[rw]*gns[rsize*n + rw]*summed[2*(lMax+1)*l*rsize + 2*m*rsize + 2*rw    ]; // Re
          Cs[2*(lMax+1)*(lMax+1)*n + l*2*(lMax+1) + 2*m + 1] += rw2[rw]*ws[rw]*gns[rsize*n + rw]*summed[2*(lMax+1)*l*rsize + 2*m*rsize + 2*rw + 1]; //Im

        }
     }
      }
    }
  isCenter[0] = 0;
}
//=========================================================
void accumC(double* Cts, double* Cs, int lMax, int gnsize, int typeI){

    for(int n = 0; n < gnsize; n++){
      for(int l = 0; l < lMax+1; l++){
        for(int m = 0; m < l+1; m++){//l+1

          Cts[2*typeI*(lMax+1)*(lMax+1)*gnsize +2*(lMax+1)*(lMax+1)*n + l*2*(lMax+1) + 2*m    ] = Cs[2*(lMax+1)*(lMax+1)*n + l*2*(lMax+1) + 2*m    ];
          Cts[2*typeI*(lMax+1)*(lMax+1)*gnsize +2*(lMax+1)*(lMax+1)*n + l*2*(lMax+1) + 2*m + 1] = Cs[2*(lMax+1)*(lMax+1)*n + l*2*(lMax+1) + 2*m + 1];

        }
      }
    }
}
//=========================================================
void getPs(double* Ps, double* Cts,  int Nt, int lMax, int gnsize){
  int NN = ((gnsize+1)*gnsize)/2;  int TT = ((Nt+1)*Nt)/2;
  int nshift = 0;
  for(int i = 0; i <TT*(lMax+1)*NN; i++){Ps[i] = 0.0;}
  int tshift = 0;

  for(int t1 = 0; t1 < Nt; t1++){
    for(int t2 = t1; t2 < Nt; t2++){
      for(int l = 0; l < lMax+1; l++){

      nshift = 0;
      for(int n = 0; n < gnsize; n++){
        for(int nd = n; nd < gnsize; nd++){
           for(int m = 0; m < l+1; m++){//l+1
              if(m==0){
                Ps[tshift*(lMax+1)*NN + l*NN + nshift ]
                 +=  Cts[2*t1*(lMax+1)*(lMax+1)*gnsize + 2*(lMax+1)*(lMax+1)*n  + l*2*(lMax+1)] // m=0
                     *Cts[2*t2*(lMax+1)*(lMax+1)*gnsize + 2*(lMax+1)*(lMax+1)*nd + l*2*(lMax+1)]; // m=0
              }else{

                Ps[tshift*(lMax+1)*NN + l*NN + nshift]
                 +=  2*(Cts[2*t1*(lMax+1)*(lMax+1)*gnsize + 2*(lMax+1)*(lMax+1)*n  + l*2*(lMax+1) + 2*m]
                      *Cts[2*t2*(lMax+1)*(lMax+1)*gnsize + 2*(lMax+1)*(lMax+1)*nd + l*2*(lMax+1) + 2*m]
		   + Cts[2*t1*(lMax+1)*(lMax+1)*gnsize + 2*(lMax+1)*(lMax+1)*n  + l*2*(lMax+1) + 2*m + 1]
                      *Cts[2*t2*(lMax+1)*(lMax+1)*gnsize + 2*(lMax+1)*(lMax+1)*nd + l*2*(lMax+1) + 2*m + 1]);
              }
            }

            nshift++;
          }
        }
      }
              tshift++;
    }
  }

}
//=========================================================
void accumP(double* Phs, double* Ps, int Nt, int lMax, int gnsize, double rCut2, int Ihpos){
  int tshift=0;
  int NN = ((gnsize+1)*gnsize)/2;
  int TT = ((Nt+1)*Nt)/2;
  for(int t1 = 0; t1 < Nt; t1++){
    for(int t2 = t1; t2 < Nt; t2++){
      for(int l = 0; l < lMax+1; l++){
       int nshift=0;

        // The power spectrum is multiplied by an l-dependent prefactor that comes
        // from the normalization of the Wigner D matrices. This prefactor is
        // mentioned in the errata of the original SOAP paper: On representing
        // chemical environments, Phys. Rev. B 87, 184115 (2013). Here the square
        // root of the prefactor in the dot-product kernel is used, so that after a
        // possible dot-product the full prefactor is recovered.
        double prefactor = PI*sqrt(8.0/(2.0*l+1.0));

        for(int n = 0; n < gnsize; n++){
          for(int nd = n; nd < gnsize; nd++){
            Phs[Ihpos*TT*(lMax+1)*NN + tshift*(lMax+1)*NN + l*NN + nshift] = prefactor*39.478417604*rCut2*Ps[tshift*(lMax+1)*NN + l*NN + nshift];// 16*9.869604401089358*Ps[tshift*(lMax+1)*NN + l*NN + nshift];
            nshift++;
          }
        }
      }
      tshift++;
    }
  }
}
//=========================================================
//=========================================================
//=========================================================
double* soap(double* c, double* Apos,double* Hpos, int* typeNs, double rCut, int totalAN,int Nt,int gnsize, int lMax, int Hs, double alpha, double* rw, double* gss);
double* soap(double* c, double* Apos,double* Hpos, int* typeNs, double rCut, int totalAN,int Nt,int gnsize, int lMax, int Hs, double alpha, double* rw, double* gss){
// everything same except last three

  double* cf = factorListSet();
  int* isCenter = malloc( sizeof(int) );
  isCenter[0] = 0;

  int rsize = 100; // constant
  double rCut2 = rCut*rCut;

  double* x    = tot  double* y    = tot  double* z    = tot double* xNow    = tot double* yNow    = tot double* zNow    = tot
  double* ris  = tot double* oOri = tot

  double* ws  = getws();
  double* oOr = getoOr(rw, rsize);  double* rw2 = getrw2(rw, rsize);

  double* oO4arri = totrs  double* minExp = totrs double* pluExp = totrs

  int Asize = 0;
  double* Cs = (double*) malloc(2*sd*(lMax+1)*(lMax+1)*gnsize);
  double* Cts = (double*) malloc(2*sd*(lMax+1)*(lMax+1)*gnsize*Nt);
  double* Ps = (double*) malloc((Nt*(Nt+1))/2*sd*(lMax+1)*((gnsize+1)*gnsize)/2);
  int icount;

  for(int Ihpos = 0; Ihpos < Hs; Ihpos++){
    for(int Itype = 0; Itype < Nt; Itype++){

        double* Ylmi; double* Flir; double* summed;
        isCenter[0] = 0;

        icount = getFilteredPos(x, y, z,xNow,yNow,zNow,ris,rw,rCut, oOri, oO4arri, minExp, pluExp,isCenter, alpha, Apos, Hpos,typeNs, rsize, Ihpos, Itype);

        Flir   = getFlir(oO4arri, ris, minExp, pluExp, icount, rsize, lMax);
        Ylmi   = getYlmi(xNow, yNow, zNow, oOri,cf,icount, lMax);
        summed = getIntegrand(Flir, Ylmi, rsize, icount, lMax);

        getC(Cs, ws, rw2, gss, summed, rCut,lMax, rsize, gnsize,isCenter,alpha);
        accumC(Cts, Cs, lMax, gnsize, Itype);

        free(Flir); free(Ylmi); free(summed);

    }

    getPs(Ps, Cts,  Nt, lMax, gnsize);
    accumP(c, Ps, Nt, lMax, gnsize,rCut2, Ihpos);
  }

  free(cf);

  free(x);  free(y);    free(z);    free(xNow);    free(yNow);    free(zNow);
  free(ris);  free(oOri);

  free(ws);
  free(oOr);  free(rw2) ;

  free(oO4arri); free(minExp); free(pluExp);
  free(Cs) ;
  free(Cts) ;
  free(Ps) ;
//  return Phs;
}
//=========================================================
//=========================================================
//=========================================================
